|
|
Event Detection - Muscular Activations (EMG) |
| Tags | detect☁emg☁tkeo |
Skeletal muscle activation is, in normal conditions, a voluntary process triggered by a nervous impulse that propagates along motor neurons until the desired muscle.
When the nervous impulse reaches sarcolemma (muscle fiber membrane) the depolarisation/repolarisation continues and the changes in membrane potential can be monitored with specialised sensors placed at skin surface.
For contracting a muscle, a large set of motor units needs to be activated, so that the acquired EMG signal is the sum of their elementary potential changes. Because of this "summation" process, EMG seems to be a little "anarchic", and the essence of EMG signal processing is in study the activation zones (bursts).
So, burst detection is an important processing step, which can be achieved by single or double threshold algorithm, generally preceded by a smoothing phase.
In this Jupyter Notebook it will be presented a single threshold algorithm, that includes the Teager-Kaiser Energy Operator (TKEO) in his implementation.
1 - Importation of the needed packages and definition of auxiliary functions
# biosignalsnotebooks python package
import biosignalsnotebooks as bsnb
# Numpy package is dedicated to simplify the work (operations between) with arrays/lists
from numpy import cumsum, concatenate, zeros, linspace, average, power, absolute, mean, std, max, array, diff, where
# Scientific packages
from scipy.signal import butter, lfilter
from scipy.stats import linregress
2 - Load of acquired EMG data
# Load of data.
data, header = bsnb.load_signal("emg_bursts", get_header=True)
3 - Identification of mac address of the device and the channel used during acquisition
channel = list(data.keys())[0]
4 - Storage of sampling rate and acquired data inside variables
# Sampling rate and acquired data
sr = header["sampling rate"]
device = header["device"]
# Signal Samples
signal = data[channel]
time = bsnb.generate_time(signal)
5 - Binarisation of EMG signal
5.1 - Preprocessing Steps
# [Baseline Removal]
pre_pro_signal = signal - average(signal)
# [Signal Filtering]
low_cutoff = 10 # Hz
high_cutoff = 300 # Hz
# Application of the signal to the filter.
pre_pro_signal = bsnb.aux_functions._butter_bandpass_filter(pre_pro_signal, low_cutoff, high_cutoff, sr)
5.2 - Application of TKEO operator
\begin{equation} TKEO[i] = \begin{cases} EMG_{original}[i], & \mbox{if } i=0 \mbox{ or } i=N-1 \\ EMG_{original}[i]^2 - (EMG_{original}[i + 1] \times EMG_{original}[i - 1]), & \mbox{otherwise}\end{cases} \end{equation} ... being $N$ the number of acquired samples.# [Application of TKEO Operator]
tkeo = []
for i in range(0, len(pre_pro_signal)):
if i == 0 or i == len(pre_pro_signal) - 1:
tkeo.append(pre_pro_signal[i])
else:
tkeo.append(power(pre_pro_signal[i], 2) - (pre_pro_signal[i + 1] * pre_pro_signal[i - 1]))
5.3 - Smoothing Phase
5.3.1 - Definition of Constants
# Smoothing level [Size of sliding window used during the moving average process (a function of sampling frequency)]
smoothing_level_perc = 20 # Percentage.
smoothing_level = int((smoothing_level_perc / 100) * sr)
5.3.2 - Signal Rectification
# [Signal Rectification]
rect_signal = absolute(tkeo)
5.3.3 - Application of the rectified signal to a first smoothing stage
# [First Moving Average Filter]
rect_signal = bsnb.aux_functions._moving_average(rect_signal, sr / 10)
5.3.4 - Application of the rectified signal to a second smoothing stage
# [Second Smoothing Phase]
smooth_signal = []
for i in range(0, len(rect_signal)):
if smoothing_level < i < len(rect_signal) - smoothing_level:
smooth_signal.append(mean(rect_signal[i - smoothing_level:i + smoothing_level]))
else:
smooth_signal.append(0)
5.4 - Definition of the detection threshold
# [Threshold]
avg_pre_pro_signal = average(pre_pro_signal)
std_pre_pro_signal = std(pre_pro_signal)
Accordingly to the method proposed by
Solnik
, threshold value can be defined as:
being $\mu_{\scriptsize EMG}$ the average EMG value, $\sigma_{\scriptsize EMG}$ his standard deviation and $h$ a variable that defines the threshold level.
To ensure that threshold level 100% is not bigger than the smooth_signal and level 0 % is not smaller than the smooth_signal we need to define a normalisation regression function.
# Regression function.
def normReg(thresholdLevel):
threshold_0_perc_level = (- avg_pre_pro_signal) / float(std_pre_pro_signal)
threshold_100_perc_level = (max(smooth_signal) - avg_pre_pro_signal) / float(std_pre_pro_signal)
m, b = linregress([0, 100], [threshold_0_perc_level, threshold_100_perc_level])[:2]
return m * thresholdLevel + b
Calculation of two threshold values
# Chosen Threshold Level (Example with two extreme values)
threshold_level = 10 # % Relative to the average value of the smoothed signal
threshold_level_norm_10 = normReg(threshold_level)
threshold_level = 80 # % Relative to the average value of the smoothed signal
threshold_level_norm_80 = normReg(threshold_level)
threshold_10 = avg_pre_pro_signal + threshold_level_norm_10 * std_pre_pro_signal
threshold_80 = avg_pre_pro_signal + threshold_level_norm_80 * std_pre_pro_signal
The threshold level of 10 % is chosen for our application, because, as can be seen in the previous figure, none activation period is completely below the threshold line.
5.5 - Binarisation of the smoothed signal
# Generation of a square wave reflecting the activation and inactivation periods.
binary_signal = []
for i in range(0, len(time)):
if smooth_signal[i] >= threshold_10:
binary_signal.append(1)
else:
binary_signal.append(0)
5.6 - Begin and end of activation periods
All upward transitions (0 to 1) define the beginning of an activation period and all downward transitions (1 to 0) establishes his enddiff_signal = diff(binary_signal)
act_begin = where(diff_signal == 1)[0]
act_end = where(diff_signal == -1)[0]
This procedure can be automatically done by detect_emg_activations function in detect module of biosignalsnotebooks package
activation_data = bsnb.detect_emg_activations(signal, sr, smooth_level=20, threshold_level=10, time_units=True, volts=False, resolution=None, device=device, plot_result=True)
As described on the intro, electromyographic (EMG) signals are generated through voluntary actions of the subject, in contrast with electrocardiographic signals.
So, due to the voluntary nature, EMG signal is not being formed uninterruptedly and between muscular activations there are inactivation periods, consisting mostly in noise, which we want to avoid during our analysis.
With the steps described on the current Jupyter Notebook , user will be in possession of an important tool to start his EMG analysis.
We hope that you have enjoyed this guide.
biosignalsnotebooks
is an environment in continuous expansion, so don"t stop your journey and learn more with the remaining
Notebooks
!